#### FIGURE 4 ####
library(tidyverse)
library(xtable)

setwd("data/analysis/")

# function to extract ATEs
Get_ATE <- function(data, subsample = c(0,1)) {
      out <- tidy(lm_robust(formula(data$FORMULA[1]), data, 
                            subset = TNRFU %in% subsample))
      out %>% slice_tail() %>%
            mutate(StudyID = data$StudyId[1])
}

# import results from heterogeneity tests
d <- read.csv("heterogeneity_stats.csv")

# Create figure
ggplot(d, aes(x= reorder(Study, r.tau.NRFUhat, sum), y=r.tau.NRFUhat)) + 
      geom_pointrange(aes(shape = factor(d.ATE.sig), 
                          ymin = r.tau.NRFUhat.LB, 
                          ymax = r.tau.NRFUhat.UB), 
                      size = .5) + 
      geom_hline(yintercept = 0, linetype = 2) + 
      xlab("Study Number") +
      ylab("Correlation Estimate") +
      scale_shape_manual(name = "Heterogeneity Detected?",
                         values = c(20, 23)) + 
      theme_bw() +
      theme(legend.position = "top", axis.text = element_text(size = 12))

# save
ggsave(file = "../../results/Figure_4.pdf")

#### TABLE I5 ####
# table of correlations
d |> select(Study, r.tau.NRFUhat, r.tau.NRFUhat.LB, r.tau.NRFUhat.UB,
            r.tau.NRFUhat.p) |>
      mutate(Study = as.integer(Study),
             r.tau.NRFUhat = round(r.tau.NRFUhat, 3),
             r.tau.NRFUhat.LB = round(r.tau.NRFUhat.LB, 3),
             r.tau.NRFUhat.UB = round(r.tau.NRFUhat.UB, 3),
             r.tau.NRFUhat.p = round(r.tau.NRFUhat.p, 3)) |>
      xtable() |>
      print(include.rownames = FALSE) |>
      write_file(file = "../../results/Table_I5.tex")


